require(data.table)
library("dplyr")
library(ggplot2)
##############################
############################  Figure 3e boxplot
##############################

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\...\\NFood_replication\\"

crops = c("Corn", "Wheat", "Soybean", "Cotton", "Rice", "Sugarcane", "Tropical Corn", "Tropical Soybean")


#******************************************************************************
#=============================Statistics and Figures===========================
#******************************************************************************

#================load data=================
#load dif.8crop.wm prepared by Data.Fig3-4.R
load(paste0(wd1,"/data/RData/Fig3-4.globalmean-yieldchanges.Rdata")) #including all grids with crop failures to calc modeled total change in yield under different scenarios (area>1000ha)

#load MLR results with confidence intervals from Bootstrap resampling (by Data3.MLR-log-bootstrap-resampling.R)
load(paste0(wd1,"/data/RData/MLR-pergrid-log-wtm-effect-Bootstrap-1000.Rdata"))


#=========aggregate per-crop effects from resampling to global total effects===========
eff.boots.pct <-  eff.log.boots[dif.8crop.wm,
                                .( Level, T.eff=(exp(T.log)-1)*100, R.eff=(exp(R.log)-1)*100, M.eff=(exp(M.log)-1)*100,
                                   RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log)-1)*100, P.eff=(exp(P.log)-1)*100, RH.eff=(exp(RH.log)-1)*100,
                                   co2.eff=(exp(i.co2.log)-1)*100, luc.eff=(exp(i.luc.log)-1)*100, #note: i.luc.log from dif.8crop (all based on SSP5 surface data) is still not the total luc effect because SSP2 have different crop grid cells than SSP5
                                   Tot.eff=(exp(Tot.log)-1)*100, Y.mod.dif=(exp(i.Y.mod.dif)-1)*100, 
                                   Tot.eff2=(exp(Tot.log2)-1)*100, Y.mod.dif2=(exp(i.Y.mod.dif2)-1)*100, #wtm of per-grid difference i.Y.mod.dif is higher than differene of global wtm yield (yield-yield0)/yield0*100
                                   Tot.eff3=(exp(Tot.log3)-1)*100, #per grid luc.log from dif.8crop.s included in tot.log3 (8.4%), larger than adding per crop i.luc.log from dif.8crop to tot.eff3 (7.4%)
                                   Y.mod.dif3=(exp(i.Y.mod.dif3)-1)*100, #per grid wtm modeled change from dif.8crop, slightly higher than (exp(i.Y.mod.dif+i.co2.log+i.luc.log)-1)*100
                                   yield=i.yield, yield1=i.yield1, yield0=i.yield0, area=i.area ,yield2=i.yield2, area2=i.area2 #prod=i.prod, prod0=i.prod0, prod1=i.prod1,
                                ), by=.EACHI, #on=.(Scenario,Year)]
                                on=.(Scenario,Crop,Irrig,Year)]
eff.boots.pct <- eff.boots.pct[!is.na(Level)]



#========Barplot to show average of mean effects and 95% confidence interval during 2075-2099 for each crop type

#====merge % effects of 8 crops to 6 crops
eff.8crop <- copy(eff.boots.pct) #not smooth
eff.8crop[Crop=="Tropical Corn", Crop:="Corn"]
eff.8crop[Crop=="Tropical Soybean", Crop:="Soybean"]

#merge tropical corn/soy to temperate corn/soy, weight by area*base yield
eff.6crop <- rbind(eff.8crop[Scenario=="RCP45 - RCP85",
                                .(T.eff=weighted.mean(T.eff, yield1*area, na.rm = T), R.eff=weighted.mean(R.eff, yield1*area, na.rm = T),M.eff=weighted.mean(M.eff, yield1*area, na.rm = T),
                                  RD.eff=weighted.mean(RD.eff, yield1*area, na.rm = T),P.eff=weighted.mean(P.eff, yield1*area, na.rm = T),
                                  RI.eff=weighted.mean(RI.eff, yield1*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield1*area, na.rm = T),
                                  Tot.eff=weighted.mean(Tot.eff, yield1*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield1*area, na.rm = T),
                                  Tot.eff2=weighted.mean(Tot.eff2, yield0*area, na.rm = T), Y.mod.dif2=weighted.mean(Y.mod.dif2, yield0*area, na.rm = T), #climate+co2 effect
                                  Tot.eff3=weighted.mean(Tot.eff3, yield0*area, na.rm = T), Y.mod.dif3=weighted.mean(Y.mod.dif3, yield0*area, na.rm = T), #climate+co2+LUC
                                  co2.eff=weighted.mean(co2.eff, yield0*area, na.rm = T), luc.eff=weighted.mean(luc.eff, yield*area, na.rm = T), #luc.log = log(i.yield2)-log(i.yield)
                                  yield0=weighted.mean(yield0, area, na.rm=T),yield1=weighted.mean(yield1, area, na.rm=T),
                                  yield=weighted.mean(yield, area, na.rm=T), yield2=weighted.mean(yield2, area2, na.rm=T), area2=sum(area2, na.rm = T),
                                  area=sum(area, na.rm = T)), by =.(Scenario,Crop,Irrig,Level,Year)],
                       eff.8crop[Scenario!="RCP45 - RCP85",
                                .(T.eff=weighted.mean(T.eff, yield0*area, na.rm = T), R.eff=weighted.mean(R.eff, yield0*area, na.rm = T),M.eff=weighted.mean(M.eff, yield0*area, na.rm = T),
                                  RD.eff=weighted.mean(RD.eff, yield0*area, na.rm = T),P.eff=weighted.mean(P.eff, yield0*area, na.rm = T),
                                  RI.eff=weighted.mean(RI.eff, yield0*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield0*area, na.rm = T),
                                  Tot.eff=weighted.mean(Tot.eff, yield0*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield0*area),
                                  yield=weighted.mean(yield, area, na.rm=T), yield0=weighted.mean(yield0, area, na.rm=T),
                                  area=sum(area, na.rm = T)),  by =.(Scenario,Crop,Irrig,Level,Year)], fill=TRUE)
eff.6crop[Scenario!="RCP45 - RCP85", co2.eff:=NA]
eff.6crop[Scenario!="RCP45 - RCP85", luc.eff:=NA]

eff.6crop.25yr <- eff.6crop[Year >=2075 & Year<=2099, lapply(.SD, mean, na.rm=T), by =.(Scenario,Crop,Irrig,Level)]
unique(eff.6crop.25yr$Crop)


#========option 1 (final): show only climate+CO2 effects for ER
eff.6crop3 <- copy(eff.6crop.25yr)
eff.6crop3 <- eff.6crop3[Scenario=="RCP45 - RCP85", Tot.eff:=Tot.eff2] #including co2.eff for ER
eff.6crop3 <- eff.6crop3[Scenario=="RCP45 - RCP85", Y.mod.dif:=Y.mod.dif2] #including co2.eff for ER

#average effects per crop at the end of 21st century (2087 is the mean of 2075-2099)
eff.6crop3.t <- eff.6crop3[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Scenario,Crop,Irrig)]
names(eff.6crop3.t)[4:7] <- c("Variable", "Mean", "Conf.L",  "Conf.U")

eff.6crop3[ Level=="Mean", Tot.eff, by=.(Scenario,Crop,Irrig)]
eff.6crop3[ Level=="Mean" & Irrig==FALSE & Scenario=="SAI - RCP85", .(RH.eff,P.eff), by=.(Scenario,Crop,Irrig)]

#report the crop-specific error term in Figure 4 (average cross scenarios)
eff.6crop3[Level=="Mean", .(Residual=mean(Tot.eff-Y.mod.dif)), by=.(Crop,Irrig)] #difference of statistical estimate and modeled change
#residual errors: -1.0%, -3%, -1.8%, 0.6%, 1.0%, 4.0% for rainfed maize, sugarcane, wheat, rice, soy, and cotton, and
#-1.5%, 0.1%, 0.1%, -0.04%, -0.4%, 0.1% for irrigated maize, sugarcane, wheat, rice, soy, and cotton, respectively




#===Final style: include CO2 effects for ER
library(RColorBrewer)
library("ggsci")
cbPalette <- pal_npg("nrc", alpha = 1)(10) #nature journal colors
# cbPalette <- pal_aaas()(10) #AAAS journal colors
# cbPalette <- pal_lancet()(9) #lancet colors
barplot(height=rep(1,length(cbPalette)),col=cbPalette) #display colors
cropcolor <- brewer.pal(n =6, "Set2") 
cropcolor <- brewer.pal(n =6, "Dark2")

#======final color for crops=====
library(viridis)
cropcolor <- viridis(18)[seq(7,18,2)]
barplot(height=rep(1,length(cropcolor)),col=cropcolor) 


#==rainfed crop
levs2 <- c("T.eff",  "RD.eff", "RI.eff", "P.eff", "RH.eff", "co2.eff", "Tot.eff", "Y.mod.dif")#
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('P'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('RH'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(CLM5)')), ' ')))

eff.6crop3.t[Scenario=="RCP45 - RCP85", Scenario:="Emissions reduction"]
scenario <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "Emissions reduction") #for facet_grid_sc
levs1 <- c("Corn", "Sugarcane", "Wheat", "Rice", "Soybean", "Cotton")
cropname = c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")

library(facetscales)
scales_y <- list(  #set different scales with facet_grid_sc
  "SAI - RCP85" = scale_y_continuous("Change in per-crop yield (%)", limits = c(-25, 53), breaks = c(-25,0,25,50)),
  "MSB - RCP85" = scale_y_continuous(limits = c(-25, 53), breaks = c(-25,0,25,50)),
  "CCT - RCP85" = scale_y_continuous(limits = c(-25, 53), breaks = c(-25,0,25,50)),
  "Emissions reduction" = scale_y_continuous(limits = c(-50, 53), breaks = c(-50,-25,0,25,50)) )

df <-  eff.6crop3.t[Variable%in%levs2 & Irrig==FALSE]
p <- df %>%  #year 2094 is the 10 year moving average of 2090-2099
  mutate(Scenario = factor(Scenario, levels=scenario, labels=scenario)) %>%
  mutate(Crop = factor(Crop, levels=levs1, labels=cropname)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  ggplot(aes(x=Variable, y=Mean, fill=Crop))+ # , alpha=Irrig)) +
  geom_bar(stat="identity", width=0.75, position=position_dodge(width=0.75)) + #position=position_dodge(width=400), binwidth=500) +
  geom_errorbar(aes(ymin=Conf.L, ymax=Conf.U), width=.2,
                position=position_dodge(width=.75)) + #smaller width increase the gap between dodge'd bin groups
  #facet_wrap(~Scenario, ncol=1, labeller= label_parsed) +
  facet_grid_sc(rows=vars(Scenario), space="free_y", scales = list(y = scales_y)) + #set different scales for facet plots
  geom_vline(aes(xintercept=6.5),linetype="dashed") + #rainfed
  scale_x_discrete(name="", labels=varname) + #label_parsed only works with facet_wrap or facet_grid
  #scale_y_continuous("Change in per-crop yield (%)", limits=c(-50, 55), breaks=c(-50,-25,0,25,50)) +
  scale_fill_manual(name="", values= cropcolor[1:6], labels=cropname) +
  #scale_alpha_manual(name="", values= c(1.,0.3)) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=0.6, keywidth=0.6))

p + theme_bw()+ #coord_flip() +theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle(expression('Rainfed')) + # ggtitle(expression('('*bold(a)*') Rainfed')) + 
  theme(plot.title = element_text(size = 10,hjust=0.5), legend.text = element_text(size=9)) +
  theme(axis.text = element_text(size=10), strip.text = element_text(size = 10), axis.title=element_text(size=10)) +
  theme(axis.text.x = element_text(angle=0,hjust=0.5,vjust=0), axis.title.x =element_blank())+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"), #legend.position = c(0.5,-0.1))
        legend.position = "bottom")

ggsave(paste0(wd1,"/figures/main/Fig4a.svg"), device ="svg", units = "in", width = 5.2, height = 6.6, dpi = 600)


#==irrigated
levs2 <- c("T.eff",  "RD.eff", "RI.eff", "co2.eff",  "Tot.eff", "Y.mod.dif")  # Irrigated
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(CLM5)')), ' ')))

scales_y <- list(  #set different scales with facet_grid_sc
  "SAI - RCP85" = scale_y_continuous(" ", limits = c(-15, 25), breaks = c(-12,0,25)),
  "MSB - RCP85" = scale_y_continuous(limits = c(-15, 25), breaks = c(-12,0,25)),
  "CCT - RCP85" = scale_y_continuous(limits = c(-15, 25), breaks = c(-12,0,25)),
  "Emissions reduction" = scale_y_continuous(limits = c(-28, 25), breaks = c(-25,0,25)) )

df <-  eff.6crop3.t[Variable%in%levs2 & Irrig==TRUE]
p <- df %>%  #year 2094 is the 10 year moving average of 2090-2099
  mutate(Scenario = factor(Scenario, levels=scenario, labels=scenario)) %>%
  mutate(Crop = factor(Crop, levels=levs1, labels=cropname)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  ggplot(aes(x=Variable, y=Mean, fill=Crop))+ # , alpha=Irrig)) +
  geom_bar(stat="identity", width=0.75, position=position_dodge(width=0.75)) + #position=position_dodge(width=400), binwidth=500) +
  geom_errorbar(aes(ymin=Conf.L, ymax=Conf.U), width=.2,
                position=position_dodge(width=.75)) + #smaller width increase the gap between dodge'd bin groups
  #facet_wrap(~Scenario, ncol=1, labeller= label_parsed) +
  facet_grid_sc(rows=vars(Scenario), space="free_y", scales = list(y = scales_y)) + #set different scales for facet plots
  geom_vline(aes(xintercept=4.5),linetype="dashed") + #irrigated
  scale_x_discrete(name="", labels=varname) + #label_parsed only works with facet_wrap or facet_grid
  #scale_y_continuous("Change in per-crop yield (%)", limits=c(-30, 25), breaks=c(-25,0,25)) +
  scale_fill_manual(name="", values= cropcolor[1:6], labels=cropname) +
  #scale_alpha_manual(name="", values= c(1.,0.3)) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=0.6, keywidth=0.6))

p + theme_bw()+ #coord_flip() +theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle('Irrigated') + #ggtitle(expression('('*bold(b)*') Irrigated')) +
  theme(plot.title = element_text(size = 10,hjust=0.5), legend.text = element_text(size=9)) +
  theme(axis.text = element_text(size=10), strip.text = element_text(size = 10), axis.title=element_text(size=10)) +
  theme(axis.text.x = element_text(angle=0,hjust=0.5,vjust=0), axis.title.x =element_blank())+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"), #legend.position = c(0.5,-0.1))
        legend.position = "bottom")

ggsave(paste0(wd1,"/figures/main/Fig4b.svg"), device ="svg", units = "in", width = 4.2, height = 6.6, dpi = 600)





#==Supplementary considering LUC effects for ER: Climate+CO2+LUC
eff.6crop2 <- copy(eff.6crop.25yr)
eff.6crop2 <- eff.6crop2[Scenario=="RCP45 - RCP85", Tot.eff:=Tot.eff3] #including co2.eff and luc.eff
eff.6crop2 <- eff.6crop2[Scenario=="RCP45 - RCP85", Y.mod.dif:=Y.mod.dif3] #including co2.eff and luc.eff

#average effects per crop at the end of 21st century
#eff.6crop2.t <- eff.6crop2[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Scenario,Crop,Irrig)]
eff.6crop2.t <- eff.6crop2[Year==2087, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Scenario,Crop,Irrig)]
names(eff.6crop2.t)[4:7] <- c("Variable", "Mean", "Conf.L",  "Conf.U")


eff.6crop2[Year==2087 & Level=="Mean", mean(Tot.eff), by=.(Scenario)] 
eff.6crop2[Year==2087 & Level=="Mean", mean(co2.eff), by=.(Scenario)]
eff.6crop2[Year==2087 & Level=="Mean", mean(luc.eff), by=.(Scenario)]
eff.6crop2[Year==2087 & Level=="Mean", mean(Y.mod.dif), by=.(Scenario)] 
#eff.6crop2[Year >=2075 & Year<=2099 & Level=="Mean", mean(Y.mod.dif3), by=.(Scenario)]


#=== include CO2 and LUC effects

#==rainfed crop
levs2 <- c("T.eff",  "RD.eff", "RI.eff", "P.eff", "RH.eff", "co2.eff", "luc.eff", "Tot.eff", "Y.mod.dif")#
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('P'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('RH'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('LUC'), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(CLM5)')), ' ')))

scenario <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85") #for facet_grid_sc
levs3 <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85")
levs1 <- c("Corn", "Sugarcane", "Wheat", "Rice", "Soybean", "Cotton")
cropname = c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")

scales_y <- list(  #set different scales with facet_grid_sc
  "SAI - RCP85" = scale_y_continuous("Change in per-crop yield (%)", limits = c(-25, 53), breaks = c(-25,0,25,50)),
  "MSB - RCP85" = scale_y_continuous(limits = c(-25, 53), breaks = c(-25,0,25,50)),
  "CCT - RCP85" = scale_y_continuous(limits = c(-25, 53), breaks = c(-25,0,25,50)),
  "RCP45 - RCP85" = scale_y_continuous(limits = c(-50, 70), breaks = c(-50,-25,0,25,50)) )

df <-  eff.6crop2.t[Variable%in%levs2 & Irrig==FALSE]
p <- df[Scenario%in%levs3] %>%  #year 2094 is the 10 year moving average of 2090-2099
  mutate(Scenario = factor(Scenario, levels=scenario, labels=scenario)) %>%
  mutate(Crop = factor(Crop, levels=levs1, labels=cropname)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  ggplot(aes(x=Variable, y=Mean, fill=Crop))+ # , alpha=Irrig)) +
  geom_bar(stat="identity", width=0.75, position=position_dodge(width=0.75)) + #position=position_dodge(width=400), binwidth=500) +
  geom_errorbar(aes(ymin=Conf.L, ymax=Conf.U), width=.2,
                position=position_dodge(width=.75)) + #smaller width increase the gap between dodge'd bin groups
  #facet_wrap(~Scenario, ncol=1, labeller= label_parsed) +
  facet_grid_sc(rows=vars(Scenario), space="free_y", scales = list(y = scales_y)) + #set different scales for facet plots
  geom_vline(aes(xintercept=7.5),linetype="dashed") + #rainfed
  scale_x_discrete(name=NULL, labels=varname) +#label_parsed only works with facet_wrap or facet_grid
  # scale_y_continuous("Change in per-crop yield (%)")+#, limits=c(-50, 70), breaks=c(-50,-25,0,25,50)) +
  scale_fill_manual(name=NULL, values= cropcolor, labels=cropname) +
  #scale_alpha_manual(name="", values= c(1.,0.3)) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=0.6, keywidth=0.4))

p + theme_bw()+ #coord_flip() +theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle(expression('('*bold(b)*') Rainfed')) + #"(b) Partial and total climate effects on rainfed crops in 2090-2099"
  theme(plot.title = element_text(size = 13,hjust=0.5), legend.text = element_text(size=11)) +
  theme(axis.text = element_text(size=12), strip.text = element_text(size = 11), axis.title=element_text(size=13)) +
  theme(axis.text.x = element_text(angle=0,hjust=0.5,vjust=0), axis.title.x =element_blank())+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"), legend.position = "bottom")
        #legend.position = c(0.39,0.3), legend.background = element_blank())
        
ggsave(paste0(wd1,"/figures/supplementary/Fig.S3b.svg"), device ="svg", units = "in", width = 5.8, height = 6.6, dpi = 600)



#==irrigated
levs2 <- c("T.eff",  "RD.eff", "RI.eff", "co2.eff", "luc.eff", "Tot.eff", "Y.mod.dif")  # Irrigated
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('LUC'), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle('Total'), textstyle('(CLM5)')), ' ')))

# scales_y <- list(  #set different scales with facet_grid_sc
#   "SAI - RCP85" = scale_y_continuous(" ", limits = c(-15, 25), breaks = c(-12,0,25)),
#   "MSB - RCP85" = scale_y_continuous(limits = c(-15, 25), breaks = c(-12,0,25)),
#   "CCT - RCP85" = scale_y_continuous(limits = c(-15, 25), breaks = c(-12,0,25)),
#   "RCP45 - RCP85" = scale_y_continuous(limits = c(-28, 70), breaks = c(-25,0,25,50)) )

df <-  eff.6crop2.t[Variable%in%levs2 & Irrig==TRUE]
p <- df[Scenario%in%levs3] %>%  
  mutate(Scenario = factor(Scenario, levels=scenario, labels=scenario)) %>%
  mutate(Crop = factor(Crop, levels=levs1, labels=cropname)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  ggplot(aes(x=Variable, y=Mean, fill=Crop))+ # , alpha=Irrig)) +
  geom_bar(stat="identity", width=0.75, position=position_dodge(width=0.75)) + #position=position_dodge(width=400), binwidth=500) +
  geom_errorbar(aes(ymin=Conf.L, ymax=Conf.U), width=.2,
                position=position_dodge(width=.75)) + #smaller width increase the gap between dodge'd bin groups
  #facet_wrap(~Scenario, ncol=1, labeller= label_parsed) +
  facet_grid_sc(rows=vars(Scenario), space="free_y", scales = list(y = scales_y)) + #set different scales for facet plots
  geom_vline(aes(xintercept=5.5),linetype="dashed") + #irrigated
  scale_x_discrete(name="", labels=varname) + #label_parsed only works with facet_wrap or facet_grid
  # scale_y_continuous("Change in per-crop yield (%)")+#, limits=c(-50, 70), breaks=c(-50,-25,0,25,50)) +
  scale_fill_manual(name="", values= cropcolor, labels=cropname) +
  #scale_alpha_manual(name="", values= c(1.,0.3)) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=0.6, keywidth=0.6))

p + theme_bw()+ #coord_flip() +theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle(expression('('*bold(c)*') Irrigated')) +
  theme(plot.title = element_text(size = 13,hjust=0.5), legend.text = element_text(size=12)) +
  theme(axis.text = element_text(size=12), strip.text = element_text(size = 11), axis.title=element_text(size=13)) +
  theme(axis.text.x = element_text(angle=0,hjust=0.5,vjust=0), axis.title.x =element_blank())+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"), #legend.position = c(0.5,-0.1))
        legend.position = "bottom")

ggsave(paste0(wd1,"/figures/supplementary/Fig.S3c.svg"), device ="svg", units = "in", width = 4.7, height = 6.6, dpi = 600)




 



